Numerical Study of the Spin Hall Conductance in the Luttinger Model 
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We present first numerical studies of the disorder effect on the recently proposed intrinsic spin Hall 
conductance in a three dimensional (3D) lattice Luttinger model. The results show that the spin 
Hall conductance remains finite in a wide range of disorder strength, with large fluctuations. The 
disorder-configuration-averaged spin Hall conductance monotonically decreases with the increase of 
disorder strength and vanishes before the Anderson localization takes place. The finite-size effect is 
also discussed. 
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A primary goal of spintronics is to make use of spin 
degree of freedom of electrons in the future 'electronic' 
devices 0, 0] . The spin Hall effect (SHE) may be one of 
potentially effective ways to manipulate the spin trans- 
port. An extrinsic SHE generated by impurities with 
spin-orbit (SO) coupling has been previously proposed 
[3J. By scattering electrons of different spins into differ- 
ent directions, a net spin current can be established in 
the transverse direction, accompanying the charge cur- 
rent induced by an applied electric field. But usually 
the resulting spin accumulation is very weak as it cru- 
cially depends on the impurity concentration. Recently, 
a much stronger SHE due to the intrinsic SO coupling 
in clean materials has been proposed for both the 3D p- 
doped semiconductors described by the Luttinger model 
, and the two-dimensional (2D) electron gas described 
by the Rashba model [jj. Here it has been argued that 
the 'dissipationless' spin currents can be of several orders 
of magnitude larger than in the case of the extrinsic SHE. 
A signature of spin polarization observed recently in the 
2D hole gas (2DHG) and 3D n-doped semiconductors 
might be originated from the intrinsic SHE || |j| . 

However, the effect of disorder on the intrinsic SHE 
remains a highly controversial issue so far. It has been 
argued [lCI 1 1 1| that the spin current in the 2D Rashba 
model should vanish, even in the weak disorder limit, 
after considering the vertex corrections. On the other 
hand, it is shown that the vertex correction vanishes for 
the Luttinger E3 and 2DHG || models such that the 
SHE is robust in the latter models at least when disor- 
ders are weak. Clearly disorder effect is nonperturbative 
on the spin Hall transport properties, and numerical ap- 
proaches are highly desirable in order to illustrate the 
fate of spin Hall conductance (SHC) in disordered sys- 
tems. So far there have been a series of numerical works 
dealing with the SHE in the presence of disorders in the 
2D Rashba model |H El El El. These calculations 
suggest that the SHC survives finite length scales in dis- 
ordered systems with indications of its vanishing in the 
thermodynamic limit. To our knowledge no numerical 
work has been done in the 3D Luttinger model with re- 
gard to the disorder effect. 



of the SHC in the lattice Luttinger model with including 
an on-site random potential based on the Kubo formula. 
We find that the SHC at weak disorder is intrinsically 
fluctuating, similar to the quantum Hall state around 
the critical point. The distribution of the SHC over dif- 
ferent disorder configurations shows a strong symmetric 
peak with the averaged SHC located at the peak posi- 
tion. The averaged SHC remains finite in a wide range 
of disorder strengths covering the main regime of the 
metallic phase while the finite-size scaling analysis sug- 
gests that SHC can survive at larger length scales. The 
calculated SHC decreases monotonically with increasing 
disorder strength, and disappears not far before the 3D 
Anderson localization takes place. 

We start with the tight-binding version of the 3D Lut- 
tinger Hamiltonian, which can be derived from the con- 
tinuum version Q with using the replacement k v — > 
sin&v and k^ — » 2(1 — cosfc„). After a discrete Fourier 
transformation, the resulting Hamiltonian reads 



H = -J2(4cj+h.c.) + V L J2(4slc l+ „ + h.c.) 
(ij) i-y 



(1) 



where the electron annihilation operator Cj has four 
components characterized by the 'spin' index S z — 
|, \, — A, — §, respectively, and i + v {y = x,y, z) denote 
the nearest-neighbors of site i, and i + v + /i, etc., for the 
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next nearest-neighbor sites. Here Vt, — 

° ^ 71 + ^72 

sents the strength of the Luttinger spin-orbital coupling. 
We choose Vl = 0.364 as the ratio 71/72 is around 3 
in typical semiconductors |l7j . The last term accounts 
for on-site nonmangetic disorder with randomly dis- 
tributed within [-W/2,W/2]. Note that the Luttinger 
model is only a valid description of real semiconductors 
around the T points at k v — > 0, which corresponds to 
choosing the Fermi energy near the band edge in the 
present tight-binding version. 
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FIG. 1: asH vs the Fermi energy Ef in the pure system. The 
dashed curve is for 8x8x8 lattice and the solid one is for 
50 x 50 x 50 lattice with PBC. The open circles for 8 x 8 x 8 
lattice are obtained by averaging over 200 configurations of 
different BCs, which coincide with the solid curve very well. 
The inset shows P(ag H ), defined as the distribution of the 
spin Hall conductance, at 8 x 8 x 8 under different BCs with 
Ef's within the range indicated by the arrow in the main 
panel. 



lated by the Kubo formulaflj) 
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E n <E f <E„ 
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(2) 



in which Nl is the number of lattice sites, Ef denotes 
the Fermi energy, E m n is the eigen-energy, the charge 
current operator j = ev and the spin current j^ spm = 
\{v^, Su}. Here the velocity operator v as the conjugate 
operator of the position operator R = J2ia r i n i<? ( n i<r 
is the number operator at site i with spin index a), is 
defined by the standard relation v = ^ [H, R] . In the 
presence of random disorder W ^ 0, the SHC is obtained 
by averaging <J l SH over all disorder configurations, i.e., 



= (4h)- 



(3) 



In a finite-size calculation, a proper boundary condi- 
tion (BC) is necessary for diagonalizing the Hamilto- 
nian. A general (twisted) BC, e.g., ip(x + L x ,y,z) = 
e 27 ™^ x ^(;E, y, z), etc., where L x is the sample size along 
the x-direction, and <px is defined within [0, 1] with <f> x = 
corresponding to the periodic BC (PBC) along this di- 
rection. In the thermodynamic limit, a physical quantity 
should not depend on BCs. In a finite-size calculation, 
the BC averaging can be very effective in smoothening 
out finite-size fluctuations in a$H m & spin-orbit coupling 
system |16| . In principle, this procedure is not necessar- 
ily the unique one for a finite system (as one can also 
use the fixed BCs in the calculation), but smoother data 
obtained this way can allow one to make a finite-size scal- 
ing analysis and to meaningfully extrapolate the results 



consider the PBC in the pure system with e,; = 0. The 
calculated asn for a 8 x 8 x 8 lattice is shown in Fig. ^by 
the dashed curve, which quickly fluctuates, as a function 
of Ef, with finite steps due to the finite-size effect. Such 
a finite-size effect disappears when the sample size is in- 
creased to 50 x 50 x 50 (this size can only be reached for 
the pure system, where the momentum is a good quan- 
tum number, in our calculation) with the same PBC, as 
indicated by the smooth solid curve in Fig. On the 
other hand, if one averages over different BCs (over 200 
configurations) in Eq. JSJ) for the 8x8x8 lattice, the 
steps in the dashed curve can also become smoothened 
out as represented by the open circles which coincide very 
well with the solid curve obtained for the bigger lattice 
of 50 x 50 x 50 in Fig. □ 

The fluctuations in cr\ H become very large in the pres- 
ence of disorders, typically in a range of 5-10 times larger 
than the averaged value. To quantitatively describe such 
fluctuations, we shall introduce the so-called distribution 
of the SHC (DSHC), P{(T X SH ), which determines the av- 
eraged SHC, <t S h, by 



a SH = J da 1 SH P{a 1 SH )a 1 SH . 



(4) 



First, for a given Ef, we can calculate <j\ h at different 
disorder and BC configurations within a small Fermi en- 
ergy interval, say, [—2.27,-2.07] around Ef = —2.17 as 
illustrated in Fig. ^ by the arrow [here the change in 
asH(Ef) is presumably weak as a function of Ef]. Sup- 
pose that the total number of computed <Jg H 's is N in 
this range, and the number of a^ H 's at <jg H = a ± Act 
[Act = ±0.01 (e/87r)], is denoted by AN {a). Then the 



DSHC is denned as the statistic distribution of a\ H 



P(a) 



AA(er) 
N Aa ' 



(5) 



The DSHC for the pure system of a 8 x 8 x 8 lattice for 
200 different BCs is shown in the inset of Fig. ^ m which 
P(cr) is a very symmetric peak such that one may simply 
use the peak position, to determine the averaged 

<7sh instead of directly evaluating the average in Eq.(@J. 
A similar technique has been used in the quantum Hall 
effect system [lflj . 

Fig. (a) shows the DSHC at W = 3 for a 6 x 6 x 6 
lattice, with the Fermi energy Ef fixed as the same value 
as in Fig. I. Here the open squares correspond to the 
result obtained over N = 1000 configurations of random 
disorder and BCs, while the closed squares are for the 
N = 5000 configurations. Clearly the DSHC becomes 
smoother with the increase of N, whose symmetric peak 
position remains unchanged with essentially the same 
envelop. The solid curve in Fig. [21 (a) is obtained by 
averaging the DSHC at N — 1000 over a small range 
of a: [a - 8, a + 6} with 6 = 0.2 (e/8?r), defined by 
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FIG. 2: (a) The DSHC atM / = 3ona6x6x6 lattice. Open 
square denotes 1000 random disorder and BC configurations, 
and the closed square is for 5000 configurations. The solid 
curve is the averaged DSHC, P a (defined in the text), at 1000 
configurations. The inset shows P a in a larger scale, (b): P a 
at different disorder strengths. The dotted one is for W = 0; 
the solid curve: W = 0.4; the dashed curve: W = 3; the 
dash-dot curve: W = 10. Inset: the same curves in a larger 
scale. 



at N = 5000 very well and is plotted in a wider range 
of a in the inset. The calculated P a at W = 0, 0.4, 
3.0, and 10, respectively, for a 8 x 8 x 8 lattice averaged 
over 200 configurations are presented in Fig. El(b). The 
main panel focus on the neighborhood around the peaks 
of P a (<r ) and the inset illustrates the peaks in a larger 
scale, whose lineshapes are generally symmetric such that 
one may read off the value of <jsh directly from the peak 
position as discussed above. 

Now we study the sample-size dependence of the SHC 
with focusing on three disorder strengths: W = 0.4 for 
the weak disorder regime; W — 3 for the intermediate 
regime; and W — 10 for the strong disorder regime. The 
results for W = 0.4 are plotted in Fig. [5] (a) for three 
different sizes of the lattice: 6x6x6 with N = 500 (the 
solid curve), 8x8x8 with N = 200 (the dashed), and 
10 x 10 x 10 with N = 200 (the dotted). In Fig. 0(b) the 
difference between the 6x6x6 and 10 x 10 x 10 lattices 
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FIG. 3: (a). The DSHC at W = 0.4. The solid curve: a 
6x6x6 lattice with 500 configurations, the dashed curve: a 
8x8x8 lattice with 200 configurations and the dotted curve: a 
10 x 10 x 10 lattice with 200 configurations, (b). The difference 
between the DSHCs of sizes 6x6x6 and 10 x 10 x 10. (c) and 
(d) are similar to (a) and (b), for W = 3 (the right side peak) 
and W = 10 (the left side peak), respectively. The insets in 
(a) and (c) shows the DSHCs at 10 x 10 x 10 for W = 0.4 and 
10, respectively, which are fit by the function -. — , — — . 



unit for all sizes, the P a at 10 x 10 x 10 has relatively much 
longer tail such that AP a in Fig. 3(b) remains negative 
over a wide range of a\ H that is not easily seen by a 
naked eye. These results show that the peak of Pa(<Jsn) 
is significantly broadened and reduced with the increase 
of the sample size, thus the fluctuations may survive in 
the large system size limit, like near the critical point of 
quantum Hall system [lj| ■ Such large fluctuations of the 
SHC may be attributed to the nonconserved spins under 
the random scattering of disorder. However, the peak 
position remains essentially unchanged, which still well 
decides the averaged <jsh ■ On the other hand, much less 
sample-size dependence is found for W = 3 and W = 10, 
corresponding to two peaks in Fig. (c), respectively, 
where the data for different sample-sizes all coincide with 
each other. And the differences of the DSHCs between 
the 6x6x6 and 10 x 10 x 10 lattices are shown in Fig. 
El (d) which are indeed much reduced as compared the 
weak disorder case in (b). 

To further characterize the size-dependence of P a and 
csh, we use a function ^ — — — — rj— to fit P a , in which 

i S H &SH\ t C 

the typical width of the DSHC, Actsh, defined as the half 

•111 1 . If ■ / T TT T 7TT1\ /T\ _1 ll 
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FIG. 4: (a) The W dependence of osh- The solid curve with 
open squares is for 6x6x6, the dashed curve with close 
squares is for 8x8x8, and the dotted curve with open circles 
is for 10 x 10 x 10, while the dash-dot curve with closed circles 
is for 6 x 6 x 24. Inset: ctsh over a wider disorder regime, 
calculated for a 10 x 10 x 10 lattice, where W c denotes the 
critical disorder for 3D Anderson localization (see text), (b) 
The W dependence of the HWHM Aasn- The inset shows 
the details in a weak disorder regime. The notations are the 
same as in the main panel of (a). 



larger W. These results suggest that there may exist a 
characteristic length scale in the spin transport which de- 
creases with the increase of W. If this is true, then the 
extrapolation to the thermodynamic limit should be at 
least reliable in the strong disorder regime for the present 
finite-size calculation. By a simple interpolation between 
the pure case and the strong disorder case, then the SHC 
is expected to be robust over a wide range of disorder 
strength. 

In conclusion, we numerically studied the distribution 
of the spin Hall conductance and determine the SHC. 
The main result shown in Fig. 0] indicates that in the 
weak disorder regime cjsh remains almost the same as 
the value for the pure system. With the increase of 
the disorder strength, asH is reduced and terminates 
before the 3D Anderson localization takes place. Al- 
though the calculation has been performed on finite lat- 
tice sizes, through the finite-size analysis of the distribu- 
tion function of the SHC, we found that the results are 
quite size- independent, suggesting that the SHC in the 
3D Luttinger model be robust. This is in contrast to the 
vanishing behavior found in 2D electron Rashba model, 
in agreement with analytical results considering vertex 
corrections [ill 0] . 

We thank the helpful discussion with S.C. Zhang, 
M.W. Wu, and C.X. Liu. This work is partially sup- 
ported by the grants from the NSFC, ACS-PRF 41752- 
AC10, and the NSF grant/DMR-0307170. The com- 
putation of this project was performed on the HP-SC45 
Sigma-X parallel computer of ITP and ICTS, CAS. 



examples of the good fitting are shown in the insets of 
Fig. |21 (a) and (c) for a 10 x 10 x 10 lattice. In this way 
we can systematically determine both asH and the cor- 
responding AasH at different sample sizes and disorder 
strengths. The results are depicted in Figs. |H(a) and (b) 
as functions of the disorder strength W. In the weak dis- 
order regime, we see that ctsh remains almost the same as 
the pure system. With the increase of W, ctsh decreases 
monotonically and is reduced to 5% of the disorder-free 
value at W ~ 10. It becomes indistinguishable from 
zero around W ~ 14, which is quite close to the typical 
critical disorder strength, W c ~ 16 ,20], of the Ander- 
son localization in 3D systems as marked in the inset of 
Fig. |U(a). So the results suggest that the SHE always 
occurs in the delocalized regime below W c . Furthermore, 
the finite-size effect is very weak at W > 2, from 6x6x6 
to 10 x 10 x 10, with the continuous reduction of AasH 
[Fig. El (b)]. 

The overall fluctuations of the SHC and the sample- 
size dependence are the strongest in the intermediate 
regime of 0.5 < W < 2. Both effects are then monotoni- 
cally reduced at W > 2, where asn becomes weakly de- 
pendent on the sample size and the relatively small Ao~sh 
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